import pandas as pd
from prophet import Prophet
from prophet.diagnostics import cross_validation, performance_metrics
import matplotlib.pyplot as plt

# 读取Excel文件
file_path = '每小时进出港统计.xlsx'
df = pd.read_excel(file_path, engine='openpyxl')

# 确保时间戳被正确解析
df['时间'] = pd.to_datetime(df['时间'])

# 将数据转换为Prophet所需的格式 (ds, y)
df_prophet = df.rename(columns={'时间': 'ds', '进出港次数': 'y'})

# 打印数据长度信息
print(f"Total data points: {len(df_prophet)}")

# 初始化Prophet模型，并设置参数
model = Prophet(
    seasonality_mode='additive',  # 季节性模式为加法
    seasonality_prior_scale=10.0,  # 较大的季节性先验规模以适应每天的周期性
    changepoint_prior_scale=0.5,   # 中等的趋势变化点先验规模
    n_changepoints=30,             # 30个趋势变化点
    interval_width=0.8             # 70%的置信区间
)

# 添加每日季节性成分
model.add_seasonality(name='daily', period=24, fourier_order=10)

# 训练模型
model.fit(df_prophet)

# 进行交叉验证
try:
    df_cv = cross_validation(
        model,
        initial='60 days',  # 初始训练期长度为30天
        period='24 hours',  # 每隔1小时进行一次预测
        horizon='24 hours'  # 预测的时间范围为1天
    )
except ValueError as e:
    print(e)
    print("Please adjust the parameters and try again.")
else:
    # 计算性能指标
    df_p = performance_metrics(df_cv)

    # 打印完整的交叉验证结果
    print("Cross Validation Results:")
    print(df_cv)

    # 打印完整的性能指标
    print("Performance Metrics:")
    print(df_p)

    # 打印平均绝对误差
    mae = df_p['mae'].mean()
    print(f"Mean Absolute Error: {mae}")

    # 打印均方误差
    mse = df_p['mse'].mean()
    print(f"Mean Squared Error: {mse}")

    # 打印均方根误差
    rmse = df_p['rmse'].mean()
    print(f"Root Mean Squared Error: {rmse}")

    # 计算数据的基本统计量
    mean_value = df_prophet['y'].mean()
    std_value = df_prophet['y'].std()
    print(f"Mean of actual values: {mean_value}")
    print(f"Standard deviation of actual values: {std_value}")

    # 可视化真实值与预测值的对比图
    fig, ax = plt.subplots(figsize=(10, 6))
    ax.plot(df_cv['ds'], df_cv['y'], label='Actual')
    ax.plot(df_cv['ds'], df_cv['yhat'], label='Predicted', linestyle='--')
    ax.fill_between(df_cv['ds'], df_cv['yhat_lower'], df_cv['yhat_upper'], color='gray', alpha=0.2)
    ax.set_xlabel('Date')
    ax.set_ylabel('Value')
    ax.set_title('Actual vs Predicted')
    ax.legend()
    plt.show()

# 预测未来日期
future = model.make_future_dataframe(periods=24, freq='h')  # 预测未来24小时的数据
forecast = model.predict(future)

# 确保预测结果为非负数
forecast['yhat'] = forecast['yhat'].clip(lower=0)
forecast['yhat_lower'] = forecast['yhat_lower'].clip(lower=0)
forecast['yhat_upper'] = forecast['yhat_upper'].clip(lower=0)

# 打印预测结果的前几行
print("\nForecast Results:")
print(forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']])

# 将预测结果保存到CSV文件
forecast.to_csv('hourly_forecast_results.csv', index=False)

# 可视化预测结果
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(df_prophet['ds'], df_prophet['y'], label='Historical Data')
ax.plot(forecast['ds'], forecast['yhat'], label='Forecast', linestyle='--')
ax.fill_between(forecast['ds'], forecast['yhat_lower'], forecast['yhat_upper'], color='gray', alpha=0.2)
ax.set_xlabel('Date')
ax.set_ylabel('Value')
ax.set_title('Historical Data and Forecast')
ax.legend()
plt.show()
